Great circle calculations for great distance travelling¶

Introduction¶

Back in high school, the topic great distance calculation was one of my best. It was infact one of those mathematical calculations I looked forward to in tests. It was a sure hack I could do blindfolded. One particular application of great circle distance calculation is in estimating the flight path distance. Long flights are a common feature in the airline industry. But what is very uncommon, especially for a cent conscious sector as the airline industry, is using a long route where a shortcut evidently exists. Russia seemed to go against this cardinal rule, for sometime in 2022, a flight from Moscow to Basel took 9 hours in contrast to the 3 hours it would have taken had the European Union not barred its airspaces to the country's airlines.

Moscow to Basel via a curve around Europe

Think of great circle distance as the shortest distance among two points on the earth's sphere. As antagonistic as the name and explanation, great circle provides true values compared to map calculated distances which are already distorted by the flat projection of the earth. There are two methods to calculate the great circle distance: harvesine and great circle distance calculation. Both formulae can be used to calculate the earth's distance between two points --so long as the path is as straight as an arrow. Additionally, to test our formulae to the limits, we will assume the Russian plane had a brief stopover at Ankara, Turkey, so that we can apply the formula to intermediate distances.

Let's load the usual packages.

In [1]:
# Import the necessary packages
import haversine as hs
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.express as px

Get coordinates for the airports¶

To calculate the great circle distance, there has to be a starting point and an ending point. Moscow was the start point, Ankara the stopover and Basel the final destination. Great circle distance calculation works with straight paths, so we will assume the plane travelled in a similar way as the crow flies. That is, no deviation whatsoever unlike in the above map!

Here are the coordinates for Moscow's Sheremetyevo airport.

In [2]:
moscow = (55.973648, 37.412503)
moscow
Out[2]:
(55.973648, 37.412503)

Here are the coordinates for Ankara's Esenboga Airport.

In [3]:
# Get the Lat-Lon coordinates of Ankara Esenboga Airport (Ankara Esenboğa Havalimanı in Turkish)
ankara = (40.124440, 32.991672)
ankara
Out[3]:
(40.12444, 32.991672)

And finally here is their final destination, the all famous EuroAirport Basel Mulhouse Freigburg (BSL) Airport in Basel, Switzerland.

In [4]:
# Get the Lat-lon coordinates of EuroAirport Basel Mulhouse Freigburg (BSL)
basel = (47.592687, 7.530862)
basel
Out[4]:
(47.592687, 7.530862)

Calculating flight path using Haversine formula¶

Haversine formula. Hmmm. Imagine this formula was formulated in the early 19th Century. It was developed to assist seamen in making the shortest route to their destinations. Without the convenience of modern day calculators, the inventor(s) of this formula was a sure genius for it even led to the development of logarithmic tables to guide sailors. Luckily, we don't have to crank up our brains to match theirs, as Python's haversine package already does the hard bit for us. But we shall try the hard way as well, albeit scaled down to a bit.

In [5]:
from IPython import display
display.Image("E:/documents/gis800_articles/jupyter/great-circle/harversine-formula.jpg")
Out[5]:

The below code line will calculate the Moscow-Ankara distance using the haversine method.

In [6]:
# Calculate distance from Moscow to Ankara
moscowAnkaraDistance = hs.haversine(moscow, ankara, unit= 'km')
moscowAnkaraDistance
Out[6]:
1791.8111208347666

The above distance from Moscow to Turkey's administrative capital is akin to travelling from Mombasa to Kisumu twice.

Let's see the total distance from Ankara to Basel, Switzerland.

In [7]:
# Calculate distance from Ankara to Basel
ankaraBaselDistance = hs.haversine(ankara, basel, unit='km')
ankaraBaselDistance
Out[7]:
2191.225246082135

From Turkey to Switzerland is roughly the same distance from Nairobi to Mombasa five times. To find the total distance the plane travelled from Moscow to Basel via Ankara, assuming no meandering in between, we will add the results of the two haversine calculations above.

In [8]:
# Find total distance from Moscow to Basel (Moscow - Ankara - Basel)
moscowAnkaraBaselDistance = print(moscowAnkaraDistance + ankaraBaselDistance)
moscowAnkaraBaselDistance

3983.036366916902

The two flight journeys fall short of many long distance flights but whoa! that total is slightly just half of the earth's radius which is 6, 378 km.

We had mentioned earlier that logarithmic tables were made to guide seamen on using the haversine formula. We shall also devise a way to validate our findings. That is, hardcoding the entire formula, yes, the entire formula in the multicolored graphic above. If they could do it by hand, can't we also do it by code, like seriously?

Hardcoding the haversine formula¶

To make sense of the symbols, we will use a color coded image to show where each formula fits in the class gcdist we've created below. gcdist shall be the function that we shall use to calculate the distance between two coordinates, only that the entire arithmetic methodology has been hard-coded.

In [9]:
from IPython import display
display.Image("E:/documents/gis800_articles/jupyter/great-circle/harversine-formula.jpg")
Out[9]:

Below is the class function gcdist.

In [10]:
# Import the math module
# This is the function for the Harvesine formula
import math

def gcdist(lat1, lon1, lat2, lon2):
    # Radius of earth at equator
    R = 6378.137
    
    # Our formula requires we conver all degrees to radians
    lat1 = math.radians(lat1)
    lat2 = math.radians(lat2)
    lon1 = math.radians(lon1)
    lon2 = math.radians(lon2)
    
    latSpan = lat1 - lat2
    lonSpan = lon1 - lon2
    
    a = math.sin(latSpan / 2) ** 2
    b = math.cos(lat1)
    c = math.cos(lat2)
    d = math.sin(lonSpan / 2) ** 2
    
    dist = 2 * R * math.asin(math.sqrt(a + b * c * d))
    
    return dist

The color codes of the functions below represent which part of the haversine formula they deal with. Each code line has been matched a color that references the specific part of the haversine formula it deals with.

In [11]:
from IPython import display
display.Image("E:/documents/gis800_articles/jupyter/great-circle/haversine-variables.jpg")
Out[11]:

Using the very long formula above, let's calculate the distance from Moscow to Basel via Ankara. If the total distance from this is close to that derived from the earlier haversine package, thumbs up to us.

In [12]:
# Calculate Moscow-Ankara distance using the harvesine function created above
lat1 = 55.973648
lon1 = 37.412503
lat2 = 40.124440
lon2 = 32.991672
moscowAnkaraHarv = gcdist(lat1, lon1, lat2, lon2)
print(moscowAnkaraHarv)

1793.815887808489

Not much difference. But will the same ring true for Ankara to Basel travel?

In [13]:
# Calculate Ankara-Basel distance using the harvesine function created above
lat3 = 40.124440
lon3 = 32.991672
lat4 = 47.592687
lon4 = 7.530862
ankaraBaselHarv = gcdist(lat3, lon3, lat4, lon4)
print(ankaraBaselHarv)

2193.676897349533

Gracia! Let's add both calculated distances.

In [14]:
moscowAnkaraBaselHarv = moscowAnkaraHarv + ankaraBaselHarv
print(moscowAnkaraBaselHarv)

3987.4927851580223

In summary, the difference in the total distance calculated from both formulaes is negligible.

Great circle distance formula¶

The following is a different method of calculating great distances. It's almost similar to the haversine formula except it is less verbose. Follow along...

In [15]:
# Great circle distance
from IPython import display
display.Image("E:/documents/gis800_articles/jupyter/great-circle/great-circle.jpg")
Out[15]:

This formula, when called into the python environment, is calculated as shown below. We have saved the formula into a class called great_circle.

In [16]:
# Calculate Moscow Ankara distance using the Great Circle Distance formula
from math import radians, degrees, sin, cos, asin, acos, sqrt

def great_circle(lon1, lat1, lon2, lat2):
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])    
    return 6371 * (
        acos(
            sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2)
        )
    )

Below is how we used great_circle function to calculate the distance from Moscow to Ankara.

In [17]:
# Now use the Great-circle formula to calculate distance from Moscow to Ankara
lat1 = 55.973648
lon1 = 37.412503
lat2 = 40.124440
lon2 = 32.991672

moscowAnkaraGD = great_circle(lon1, lat1, lon2, lat2)
print(moscowAnkaraGD)

1791.8086458832574

The same function has been applied to find out the route from Ankara to Basel.

In [18]:
# Now use the Great-circle formula to calculate distance from Ankara to Basel
lat3 = 40.124440
lon3 = 32.991672
lat4 = 47.592687
lon4 = 7.530862

ankaraBaselGD = great_circle(lon3, lat3, lon4, lat4)
print(ankaraBaselGD)

2191.222219437098

Let's find the total.

In [19]:
moscowAnkaraBaselGDkm = moscowAnkaraGD + ankaraBaselGD
print(moscowAnkaraBaselGDkm)

3983.0308653203556

This last calculation seems like the approach the harvesine package adopted. In fact, the calculations and the results are as similar as two peas. Just a few lines were saved, and thus, in our view, this should be the go to formula if you are math fan.

Drawing a great circle map¶

Even if many cartographers were also mathematicians, think Eratosthenes, they had to produce a map somewhere. We don't like breaking tradition, especially from great minds like Eratosthenes, so we will strive to create one. Maths has to meet with map. They have always been friends anyway!

Alright. Agreed. But what kind of map? Not just any ordinary map, but an orthographic projection. Why? Our flight paths will conspicously be shown as a curved path much like how one would view the flight travel from space. We shall work with plotly in this, thanks to its simplicity. plotly visualizations works best with a database so we will create a table of our airports along with any auxillary data.

In [20]:
# Draw a great circle map
# First create a table of the airport stations and their lat lon coordinates

airports = {
    'airportName': ['sheretyevo', 'ankara', 'euroairport basel'], 
    'country': ['russia', 'turkey', 'switzerland'],
    'latitude': [55.973648, 40.124440, 47.592687],
    'longitude': [37.412503, 32.991672, 7.530862]
}

airportDf = pd.DataFrame(airports)
airportDf
Out[20]:
airportName country latitude longitude
0 sheretyevo russia 55.973648 37.412503
1 ankara turkey 40.124440 32.991672
2 euroairport basel switzerland 47.592687 7.530862

Below is the fligth path created using plotly. Feel free to experiment with other projections.

In [21]:
# Draw in Plotly the flight path
df = airportDf
fig = px.line_geo(df,
                 lat='latitude',
                 lon='longitude',
                  projection='orthographic'
                 )

fig.show()

The advantage of plotly over other plotting packages is its interactivity. If you take a peek at the top right margin, you will find other tools which allow panning, zooming and downloading your image. Not many visualization package come with tools. Back to our map, however, the flight paths seem more like straight lines rather than curved routes. If we zoom to the area of focus, the curvature of the flight paths should be more pronounced.

In [22]:
# Show flight paths zoomed into Europe
df = airportDf
fig = px.line_geo(df,
                 lat='latitude',
                 lon='longitude',
                 text='airportName',
                 projection='orthographic',
                  scope='europe',
                  title="Flight Path of Russian Plane after Europe flight ban",
                  width=1000,
                  height=600
                 )

fig.show()

That's more like it. But just for fun, we want a more pronounced curved flight path. To display a curvy path, we will extend the flight path from Basel to the other side of the globe. Think Ezeiza Airport in Argentina.

In [23]:
# Add Argentine airport - Ezeiza airport
airports2 = {
    'airportName': ['sheretyevo', 'ankara', 'euroairport basel', 'ezeiza airport'], 
    'country': ['russia', 'turkey', 'switzerland', 'argentina'],
    'latitude': [55.973648, 40.124440, 47.592687, -34.817337],
    'longitude': [37.412503, 32.991672, 7.530862, -58.531080]
}

airportDf2 = pd.DataFrame(airports2)
airportDf2
Out[23]:
airportName country latitude longitude
0 sheretyevo russia 55.973648 37.412503
1 ankara turkey 40.124440 32.991672
2 euroairport basel switzerland 47.592687 7.530862
3 ezeiza airport argentina -34.817337 -58.531080

Now an orthographic projection of this flight path should show the ellipsoidal influence of the earth on our studied flight travels.

In [24]:
df2 = airportDf2
fig = px.line_geo(df2,
                 lat='latitude',
                 lon='longitude',
                 text='airportName',
                 projection='orthographic',
                  scope='world',
                  center= {'lat': -16.015741, 'lon':  -5.734472},
                  title="Flight Path of Russian Plane after Europe flight ban",
                  #hover_data = ['country', 'latitude', 'longitude', 'airportName'],
                  hover_data = {'country', 'latitude', 'longitude', 'airportName'},
                  basemap_visible=True,
                  width=1000,
                  height=600
                 )

fig.show()

This map projection depicts the flight curve quite well. Play around with the map, try around different projections. Have fun with interactivity of plotly and the various projections!

Conclusion¶

The haversine method is a common calculation for finding out the distance between two points, assuming the object travelled in a straight line. Although it does not account for deviations, it provides an estimate of the total distance between two points. The orthographic projection on the other hand depicts the distance travelled by an object quite succinctly, in respect to the spheroid shape of the earth.

In [ ]:
In [ ]: